
set.seed(123)

rm(list = ls())

library("readr")
library("estimatr")
library("dplyr")
library("ggplot2")
library("interflex")
library("ggpubr")
library("scales")
library("texreg")
library("vtable")
library("tableone")

theme_set(theme_bw(base_family = "Times New Roman"))

########################
# Load data
########################

data <- read_csv("data.csv")

########################
# Transform data
########################

long_data <- select(data, c("ResponseId", 
                            "article1_rating", "article1_condition", 
                            "article2_rating", "article2_condition",
                            "article3_rating", "article3_condition",
                            "article4_rating", "article4_condition",
                            "article5_rating", "article5_condition",
                            "article6_rating", "article6_condition", 
                            "feeling_DPP", "feeling_KMT", "feeling_TPP", "feeling_invalid",
                            "feeling_DPP_positive", "feeling_KMT_positive", "feeling_TPP_positive",
                            "feeling_DPP_neutral", "feeling_KMT_neutral", "feeling_TPP_neutral",
                            "vote_DPP", "vote_KMT", "vote_TPP", "vote_none"))

long_data1 <- select(long_data, c("ResponseId", "article1_rating", "article1_condition", 
                                  "feeling_DPP", "feeling_KMT", "feeling_TPP", "feeling_invalid",
                                  "feeling_DPP_positive", "feeling_KMT_positive", "feeling_TPP_positive",
                                  "feeling_DPP_neutral", "feeling_KMT_neutral", "feeling_TPP_neutral",
                                  "vote_DPP", "vote_KMT", "vote_TPP", "vote_none"))
long_data2 <- select(long_data, c("ResponseId", "article2_rating", "article2_condition", 
                                  "feeling_DPP", "feeling_KMT", "feeling_TPP", "feeling_invalid",
                                  "feeling_DPP_positive", "feeling_KMT_positive", "feeling_TPP_positive",
                                  "feeling_DPP_neutral", "feeling_KMT_neutral", "feeling_TPP_neutral",
                                  "vote_DPP", "vote_KMT", "vote_TPP", "vote_none"))
long_data3 <- select(long_data, c("ResponseId", "article3_rating", "article3_condition", 
                                  "feeling_DPP", "feeling_KMT", "feeling_TPP", "feeling_invalid",
                                  "feeling_DPP_positive", "feeling_KMT_positive", "feeling_TPP_positive",
                                  "feeling_DPP_neutral", "feeling_KMT_neutral", "feeling_TPP_neutral",
                                  "vote_DPP", "vote_KMT", "vote_TPP", "vote_none"))
long_data4 <- select(long_data, c("ResponseId", "article4_rating", "article4_condition", 
                                  "feeling_DPP", "feeling_KMT", "feeling_TPP", "feeling_invalid",
                                  "feeling_DPP_positive", "feeling_KMT_positive", "feeling_TPP_positive",
                                  "feeling_DPP_neutral", "feeling_KMT_neutral", "feeling_TPP_neutral",
                                  "vote_DPP", "vote_KMT", "vote_TPP", "vote_none"))
long_data5 <- select(long_data, c("ResponseId", "article5_rating", "article5_condition", 
                                  "feeling_DPP", "feeling_KMT", "feeling_TPP", "feeling_invalid",
                                  "feeling_DPP_positive", "feeling_KMT_positive", "feeling_TPP_positive",
                                  "feeling_DPP_neutral", "feeling_KMT_neutral", "feeling_TPP_neutral",
                                  "vote_DPP", "vote_KMT", "vote_TPP", "vote_none"))
long_data6 <- select(long_data, c("ResponseId", "article6_rating", "article6_condition", 
                                  "feeling_DPP", "feeling_KMT", "feeling_TPP", "feeling_invalid",
                                  "feeling_DPP_positive", "feeling_KMT_positive", "feeling_TPP_positive",
                                  "feeling_DPP_neutral", "feeling_KMT_neutral", "feeling_TPP_neutral",
                                  "vote_DPP", "vote_KMT", "vote_TPP", "vote_none"))

long_data1 <- rename(long_data1, article_rating = article1_rating, article_condition = article1_condition)
long_data1$article <- "Article 1"

long_data2 <- rename(long_data2, article_rating = article2_rating, article_condition = article2_condition)
long_data2$article <- "Article 2"

long_data3 <- rename(long_data3, article_rating = article3_rating, article_condition = article3_condition)
long_data3$article <- "Article 3"

long_data4 <- rename(long_data4, article_rating = article4_rating, article_condition = article4_condition)
long_data4$article <- "Article 4"

long_data5 <- rename(long_data5, article_rating = article5_rating, article_condition = article5_condition)
long_data5$article <- "Article 5"

long_data6 <- rename(long_data6, article_rating = article6_rating, article_condition = article6_condition)
long_data6$article <- "Article 6"

long_data <- rbind(long_data1, long_data2, long_data3, long_data4, long_data5, long_data6)
long_data <- long_data[order(long_data$ResponseId, long_data$article),]

########################
# Summary statistics
########################

data$feeling_none <- ifelse(data$feeling_DPP_neutral == 1 & data$feeling_KMT_neutral == 1 & data$feeling_TPP_neutral == 1, 1, 0)

st(data,
   vars = c("female", "age", "education", "trust_media", "trust_socialmedia", "political_interest",
            "feeling_DPP", "feeling_KMT", "feeling_TPP",
            "vote_party", "vote_none"),
   summ = c("notNA(x)", "mean(x)", "sd(x)"),
   summ.names = c("N", "Mean", "Std. Dev."),
   labels = c("Female (0/1)", "Age (1-5)", "Education (1-7)", "Turst in media (1-4)", "Trust in social media (1-4)", "Political interest (1-4)",
              "Feeling DPP (0-10)", "Feeling KMT (0-10)", "Feeling TPP (0-10)",
              "Voted", "Did not vote"),
   numformat = formatfunc(digits = 2, nsmall = 2, big.mark = ""),
   skip.format = c("notNA(x)"),
   out = "latex", file = "results/Table_sumstat.tex"
)

CreateTableOne(vars = c("article1_rating"), strata = "article1_condition", data = data, test = FALSE)
CreateTableOne(vars = c("article2_rating"), strata = "article2_condition", data = data, test = FALSE)
CreateTableOne(vars = c("article3_rating"), strata = "article3_condition", data = data, test = FALSE)
CreateTableOne(vars = c("article4_rating"), strata = "article4_condition", data = data, test = FALSE)
CreateTableOne(vars = c("article5_rating"), strata = "article5_condition", data = data, test = FALSE)
CreateTableOne(vars = c("article6_rating"), strata = "article6_condition", data = data, test = FALSE)

CreateTableOne(vars = c("article_rating"), strata = "article_condition", data = long_data, test = FALSE)

########################
# Balance checks
########################

## Article 1

covariates <- c("female", "age", "education", "trust_media", "trust_socialmedia", "political_interest",
                "feeling_DPP", "feeling_KMT", "feeling_TPP",
                "vote_DPP", "vote_KMT", "vote_TPP", "vote_none")

balance_table <- CreateTableOne(vars = covariates, strata = "article1_condition", data = data, test = TRUE)
balance_table <- print(balance_table)

write.csv(balance_table, file = "results/Table_balance_article1.csv")

## Article 2

covariates <- c("female", "age", "education", "trust_media", "trust_socialmedia", "political_interest",
                "feeling_DPP", "feeling_KMT", "feeling_TPP",
                "vote_DPP", "vote_KMT", "vote_TPP", "vote_none")

balance_table <- CreateTableOne(vars = covariates, strata = "article2_condition", data = data)
balance_table <- print(balance_table)

write.csv(balance_table, file = "results/Table_balance_article2.csv")

## Article 3

covariates <- c("female", "age", "education", "trust_media", "trust_socialmedia", "political_interest",
                "feeling_DPP", "feeling_KMT", "feeling_TPP",
                "vote_DPP", "vote_KMT", "vote_TPP", "vote_none")

balance_table <- CreateTableOne(vars = covariates, strata = "article3_condition", data = data)
balance_table <- print(balance_table)

write.csv(balance_table, file = "results/Table_balance_article3.csv")

## Article 4

covariates <- c("female", "age", "education", "trust_media", "trust_socialmedia", "political_interest",
                "feeling_DPP", "feeling_KMT", "feeling_TPP",
                "vote_DPP", "vote_KMT", "vote_TPP", "vote_none")

balance_table <- CreateTableOne(vars = covariates, strata = "article4_condition", data = data)
balance_table <- print(balance_table)

write.csv(balance_table, file = "results/Table_balance_article4.csv")

## Article 5

covariates <- c("female", "age", "education", "trust_media", "trust_socialmedia", "political_interest",
                "feeling_DPP", "feeling_KMT", "feeling_TPP",
                "vote_DPP", "vote_KMT", "vote_TPP", "vote_none")

balance_table <- CreateTableOne(vars = covariates, strata = "article5_condition", data = data)
balance_table <- print(balance_table)

write.csv(balance_table, file = "results/Table_balance_article5.csv")

## Article 6

covariates <- c("female", "age", "education", "trust_media", "trust_socialmedia", "political_interest",
                "feeling_DPP", "feeling_KMT", "feeling_TPP",
                "vote_DPP", "vote_KMT", "vote_TPP", "vote_none")

balance_table <- CreateTableOne(vars = covariates, strata = "article6_condition", data = data)
balance_table <- print(balance_table)

write.csv(balance_table, file = "results/Table_balance_article6.csv")

########################
# Explore
########################

st(data,
   vars = c("female", "age", "education", "trust_media", "trust_socialmedia", "political_interest",
            "feeling_DPP", "feeling_KMT", "feeling_TPP",
            "vote_party", "vote_none"),
   summ = c("notNA(x)", "mean(x)", "sd(x)"),
   summ.names = c("N", "Mean", "Std. Dev."),
   labels = c("Female (0/1)", "Age (1-5)", "Education (1-7)", "Turst in media (1-4)", "Trust in social media (1-4)", "Political interest (1-4)",
              "Feeling DPP (0-10)", "Feeling KMT (0-10)", "Feeling TPP (0-10)",
              "Voted", "Did not vote"),
   numformat = formatfunc(digits = 2, nsmall = 2, big.mark = ""),
   skip.format = c("notNA(x)"),
   out = "latex", file = "results/Table_sumstat.tex"
)

########################
# Differences between partisans and moderates
########################

data$partisans[data$feeling_DPP_positive == 1 | data$feeling_KMT_positive == 1 | data$feeling_TPP_positive == 1] <- 1
data$partisans[data$feeling_invalid == 1] <- NA
data$moderates[data$feeling_DPP_neutral == 1 & data$feeling_KMT_neutral == 1 & data$feeling_TPP_neutral] <- 1

data$respondent_type[data$partisans == 1] <- "partisans"
data$respondent_type[data$moderates == 1] <- "moderates"
data$respondent_type <- factor(data$respondent_type, levels = c("partisans", "moderates"))

mod1 <- lm_robust(female ~ respondent_type, se_type = "stata", data = data)
mod2 <- lm_robust(age ~ respondent_type, se_type = "stata", data = data)
mod3 <- lm_robust(education ~ respondent_type, se_type = "stata", data = data)
mod4 <- lm_robust(trust_media ~ respondent_type, se_type = "stata", data = data)
mod5 <- lm_robust(trust_socialmedia ~ respondent_type, se_type = "stata", data = data)
mod6 <- lm_robust(political_interest ~ respondent_type, se_type = "stata", data = data)
mod7 <- lm_robust(vote_none ~ respondent_type, se_type = "stata", data = data)
mod8 <- lm_robust(log(duration) ~ respondent_type, se_type = "stata", data = data)

texreg(list(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8),
       type = "latex",
       file = "results/Table_moderates.tex",
       digits = 3,
       stars = c(0.001, 0.01, 0.05, 0.1),
       symbol = "\\sym{+}",
       include.ci = FALSE,
       omit.coef = "(Intercept)",
       custom.header = list("Female" = 1, "Age" = 2, "Education" = 3, "Trust in media" = 4, "Trust in social media" = 5, "Political interest" = 6, "Did not vote" = 7, "Survey completion time" = 8),
       custom.model.names = c("(1)", "(2)", "(3)", "(4)", "(5)", "(6)", "(7)", "(8)"),
       custom.coef.map = list("respondent_typemoderates" = "Moderates"),
       custom.note = "")

########################
# All respondents
########################

# estimate

mod1 <- lm_robust(article_rating ~ article_condition + as.factor(article), clusters = ResponseId, se_type = "stata", data = long_data)
mod2 <- lm_robust(article_rating ~ article_condition + as.factor(article) + as.factor(ResponseId), clusters = ResponseId, se_type = "stata", data = long_data)

# table

texreg(list(mod1, mod2),
       type = "latex",
       file = "results/Table_bias_perception_all.tex",
       digits = 3,
       stars = c(0.001, 0.01, 0.05, 0.1),
       symbol = "\\sym{+}",
       include.ci = FALSE,
       omit.coef = "(Intercept)",
       custom.header = list("Article Rating" = 1:2),
       custom.model.names = c("(1)", "(2)"),
       custom.coef.map = list("article_conditionCNA" = "Central News Agency", "article_conditionCT" = "China Times", "article_conditionLT" = "Liberty Times", "article_conditionUD" = "United Daily"),
       custom.gof.rows = list("Article FEs" = c("Yes", "Yes"), "Individual FEs" = c("No", "Yes")),
       custom.note = "")

# figure

tidy_mod1 <- tidy(mod1, conf.int = TRUE) %>% mutate(model = "Article FEs")
tidy_mod2 <- tidy(mod2, conf.int = TRUE) %>% mutate(model = "Article FEs + Individual FEs")

combined_models <- rbind(tidy_mod1, tidy_mod2)
combined_models <- subset(combined_models, combined_models$term != "(Intercept)")
combined_models <- subset(combined_models, !grepl("as.factor", term))
combined_models$term <- factor(combined_models$term, levels = c("article_conditionCNA", "article_conditionCT", "article_conditionLT", "article_conditionUD"), labels = c("Central News Agency", "China Times", "Liberty Times", "United Daily"))
combined_models$conf.low.90 <- combined_models$estimate - (qnorm(0.95)*combined_models$std.error)
combined_models$conf.high.90 <- combined_models$estimate + (qnorm(0.95)*combined_models$std.error)
combined_models$model <- factor(combined_models$model, levels = c("Article FEs", "Article FEs + Individual FEs"))
combined_models$effect <- "effect"

Figure <- ggplot(combined_models, aes(x = estimate, y = effect, color = model, shape = model)) + 
  facet_wrap("term", nrow = 1) + 
  geom_point(position = position_dodge(width = 0.5), size = 2) + 
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0, position = position_dodge(width = 0.5)) +
  geom_errorbarh(aes(xmin = conf.low.90, xmax = conf.high.90), height = 0, position = position_dodge(width = 0.5), size = 1) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75, alpha = 1.0) +
  labs(x = "Estimated Treatment Effects", y = NULL, color = "Model", shape = "Model", title = NULL) +
  theme(legend.position = "bottom", axis.title.x = element_text(margin = margin(t = 5)), axis.text.x = element_blank(), axis.ticks.x = element_blank()) + 
  scale_x_continuous(limits = c(-1, 1), breaks = seq(-1, 1, length.out = 5), labels = label_number(accuracy = 0.01)) + 
  scale_color_manual(values = c("Article FEs" = "darkred", "Article FEs + Individual FEs" = "darkblue")) + 
  scale_shape_manual(values = c("Article FEs" = 16, "Article FEs + Individual FEs" = 17)) + 
  guides(shape = guide_legend(nrow = 1), color = guide_legend(nrow = 1)) + 
  coord_flip()
Figure

ggsave("results/Figure_bias_perception_all.pdf", plot = Figure, width = 10, height = 4.5, units = "in", device = cairo_pdf)

########################
# All respondents (by article)
########################

# estimate

mod1 <- lm_robust(article1_rating ~ article1_condition, se_type = "stata", data = data)
mod2 <- lm_robust(article2_rating ~ article2_condition, se_type = "stata", data = data)
mod3 <- lm_robust(article3_rating ~ article3_condition, se_type = "stata", data = data)
mod4 <- lm_robust(article4_rating ~ article4_condition, se_type = "stata", data = data)
mod5 <- lm_robust(article5_rating ~ article5_condition, se_type = "stata", data = data)
mod6 <- lm_robust(article6_rating ~ article6_condition, se_type = "stata", data = data)

# table

texreg(list(mod1, mod2, mod3, mod4, mod5, mod6),
       type = "latex",
       file = "results/Table_bias_perception_byarticle.tex",
       digits = 3,
       stars = c(0.001, 0.01, 0.05, 0.1),
       symbol = "\\sym{+}",
       include.ci = FALSE,
       omit.coef = "(Intercept)",
       custom.header = list("Article Rating" = 1:6),
       custom.model.names = c("(1)", "(2)", "(3)", "(4)", "(5)", "(6)"),
       custom.coef.map = list("article1_conditionCNA" = "Central News Agency", "article2_conditionCNA" = "Central News Agency", "article3_conditionCNA" = "Central News Agency", "article4_conditionCNA" = "Central News Agency", "article5_conditionCNA" = "Central News Agency", "article6_conditionCNA" = "Central News Agency",
                              "article1_conditionCT" = "China Times", "article2_conditionCT" = "China Times", "article3_conditionCT" = "China Times", "article4_conditionCT" = "China Times", "article5_conditionCT" = "China Times", "article6_conditionCT" = "China Times",
                              "article1_conditionLT" = "Liberty Times", "article2_conditionLT" = "Liberty Times", "article3_conditionLT" = "Liberty Times", "article4_conditionLT" = "Liberty Times", "article5_conditionLT" = "Liberty Times", "article6_conditionLT" = "Liberty Times",
                              "article1_conditionUD" = "United Daily", "article2_conditionUD" = "United Daily", "article3_conditionUD" = "United Daily", "article4_conditionUD" = "United Daily", "article5_conditionUD" = "United Daily", "article6_conditionUD" = "United Daily"),
       custom.gof.rows = list("Article" = c("1", "2", "3", "4", "5", "6")),
       custom.note = "")

# figure

tidy_mod1 <- tidy(mod1, conf.int = TRUE) %>% mutate(article = "Article 1")
tidy_mod2 <- tidy(mod2, conf.int = TRUE) %>% mutate(article = "Article 2")
tidy_mod3 <- tidy(mod3, conf.int = TRUE) %>% mutate(article = "Article 3")
tidy_mod4 <- tidy(mod4, conf.int = TRUE) %>% mutate(article = "Article 4")
tidy_mod5 <- tidy(mod5, conf.int = TRUE) %>% mutate(article = "Article 5")
tidy_mod6 <- tidy(mod6, conf.int = TRUE) %>% mutate(article = "Article 6")

combined_models <- rbind(tidy_mod1, tidy_mod2, tidy_mod3, tidy_mod4, tidy_mod5, tidy_mod6)
combined_models <- subset(combined_models, combined_models$term != "(Intercept)")
combined_models$term[grepl("conditionCNA", combined_models$term)] <- "article_conditionCNA"
combined_models$term[grepl("conditionCT", combined_models$term)] <- "article_conditionCT"
combined_models$term[grepl("conditionLT", combined_models$term)] <- "article_conditionLT"
combined_models$term[grepl("conditionUD", combined_models$term)] <- "article_conditionUD"
combined_models$term <- factor(combined_models$term, levels = c("article_conditionCNA", "article_conditionCT", "article_conditionLT", "article_conditionUD"), labels = c("Central News Agency", "China Times", "Liberty Times", "United Daily"))
combined_models$conf.low.90 <- combined_models$estimate - (qnorm(0.95)*combined_models$std.error)
combined_models$conf.high.90 <- combined_models$estimate + (qnorm(0.95)*combined_models$std.error)
combined_models$effect <- "effect"

Figure <- ggplot(combined_models, aes(x = estimate, y = effect, shape = article)) + 
  facet_wrap("term", nrow = 1) + 
  geom_point(position = position_dodge(width = 0.5), size = 2) + 
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0, position = position_dodge(width = 0.5)) +
  geom_errorbarh(aes(xmin = conf.low.90, xmax = conf.high.90), height = 0, position = position_dodge(width = 0.5), size = 1) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75, alpha = 1.0) +
  labs(x = "Estimated Treatment Effects", y = NULL, shape = "Article", title = NULL) +
  theme(legend.position = "bottom", axis.title.x = element_text(margin = margin(t = 5)), axis.text.x = element_blank(), axis.ticks.x = element_blank()) + 
  scale_x_continuous(limits = c(-1.5, 1.5), breaks = seq(-1.5, 1.5, length.out = 5), labels = label_number(accuracy = 0.01)) + 
  guides(shape = guide_legend(nrow = 1)) + 
  coord_flip()
Figure 

ggsave("results/Figure_bias_perception_byarticle.pdf", plot = Figure, width = 10, height = 4.5, units = "in", device = cairo_pdf)

########################
# Respondents by partisanship (feeling thermometer)
########################

# estimate

mod1 <- lm_robust(article_rating ~ article_condition + as.factor(article) + as.factor(ResponseId), clusters = ResponseId, se_type = "stata", data = subset(long_data, long_data$feeling_invalid != 1 & long_data$feeling_DPP_positive == 1))
mod2 <- lm_robust(article_rating ~ article_condition + as.factor(article) + as.factor(ResponseId), clusters = ResponseId, se_type = "stata", data = subset(long_data, long_data$feeling_invalid != 1 & long_data$feeling_KMT_positive == 1))
mod3 <- lm_robust(article_rating ~ article_condition + as.factor(article) + as.factor(ResponseId), clusters = ResponseId, se_type = "stata", data = subset(long_data, long_data$feeling_invalid != 1 & long_data$feeling_TPP_positive == 1))
mod4 <- lm_robust(article_rating ~ article_condition + as.factor(article) + as.factor(ResponseId), clusters = ResponseId, se_type = "stata", data = subset(long_data, long_data$feeling_DPP_neutral == 1 & long_data$feeling_KMT_neutral == 1 & long_data$feeling_TPP_neutral == 1))

# table

texreg(list(mod1, mod2, mod3, mod4),
       type = "latex",
       file = "results/Table_bias_perception_bypartisanship_feel_bothFEs.tex",
       digits = 3,
       stars = c(0.001, 0.01, 0.05, 0.1),
       symbol = "\\sym{+}",
       include.ci = FALSE,
       omit.coef = "(Intercept)",
       custom.header = list("Article Rating" = 1:4),
       custom.model.names = c("(1)", "(2)", "(3)", "(4)"),
       custom.coef.map = list("article_conditionCNA" = "Central News Agency", 
                              "article_conditionCT" = "China Times",
                              "article_conditionLT" = "Liberty Times", 
                              "article_conditionUD" = "United Daily"),
       custom.gof.rows = list("Respondents" = c("DPP supporters", "KMT supporters", "TPP supporters", "Moderates"),
                              "Article FEs" = c("Yes", "Yes", "Yes", "Yes"), "Individual FEs" = c("Yes", "Yes", "Yes", "Yes")),
       custom.note = "")

# figure

tidy_mod1 <- tidy(mod1, conf.int = TRUE) %>% mutate(respondent = "DPP supporters")
tidy_mod2 <- tidy(mod2, conf.int = TRUE) %>% mutate(respondent = "KMT supporters")
tidy_mod3 <- tidy(mod3, conf.int = TRUE) %>% mutate(respondent = "TPP supporters")
tidy_mod4 <- tidy(mod4, conf.int = TRUE) %>% mutate(respondent = "Moderates")

combined_models <- rbind(tidy_mod1, tidy_mod2, tidy_mod3, tidy_mod4)
combined_models <- subset(combined_models, combined_models$term != "(Intercept)")
combined_models <- subset(combined_models, !grepl("as.factor", term))
combined_models$term <- factor(combined_models$term, levels = c("article_conditionCNA", "article_conditionCT", "article_conditionLT", "article_conditionUD"), labels = c("Central News Agency", "China Times", "Liberty Times", "United Daily"))
combined_models$respondent <- factor(combined_models$respondent, levels = c("DPP supporters", "KMT supporters", "TPP supporters", "Moderates"))
combined_models$conf.low.90 <- combined_models$estimate - (qnorm(0.95)*combined_models$std.error)
combined_models$conf.high.90 <- combined_models$estimate + (qnorm(0.95)*combined_models$std.error)
combined_models$effect <- "effect"

Figure <- ggplot(combined_models, aes(x = estimate, y = effect, color = respondent, shape = respondent)) + 
  facet_wrap("term", nrow = 1) + 
  geom_point(position = position_dodge(width = 0.5), size = 2) + 
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0, position = position_dodge(width = 0.5)) +
  geom_errorbarh(aes(xmin = conf.low.90, xmax = conf.high.90), height = 0, position = position_dodge(width = 0.5), size = 1) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75, alpha = 1.0) +
  labs(x = "Estimated Treatment Effects", y = NULL, color = "Respondent", shape = "Respondent", title = NULL) +
  theme(legend.position = "bottom", axis.title.x = element_text(margin = margin(t = 5)), axis.text.x = element_blank(), axis.ticks.x = element_blank()) + 
  scale_x_continuous(limits = c(-2, 2), breaks = seq(-2, 2, length.out = 5), labels = label_number(accuracy = 0.01)) + 
  scale_color_manual(values = c("DPP supporters" = "darkgreen", "KMT supporters" = "darkblue", "TPP supporters" = rgb(0, 0.6, 0.6), "Moderates" = "darkgrey")) + 
  scale_shape_manual(values = c("DPP supporters" = 16, "KMT supporters" = 17, "TPP supporters" = 18, "Moderates" = 4)) + 
  guides(shape = guide_legend(nrow = 1), color = guide_legend(nrow = 1)) + 
  coord_flip()
Figure 

ggsave("results/Figure_bias_perception_bypartisanship_feel_bothFEs.pdf", plot = Figure, width = 10, height = 4.5, units = "in", device = cairo_pdf)

########################
# Respondents by partisanship (vote)
########################

# estimate

mod1 <- lm_robust(article_rating ~ article_condition + as.factor(article), clusters = ResponseId, se_type = "stata", data = subset(long_data, long_data$vote_DPP == 1))
mod2 <- lm_robust(article_rating ~ article_condition + as.factor(article), clusters = ResponseId, se_type = "stata", data = subset(long_data, long_data$vote_KMT == 1))
mod3 <- lm_robust(article_rating ~ article_condition + as.factor(article), clusters = ResponseId, se_type = "stata", data = subset(long_data, long_data$vote_TPP == 1))
mod4 <- lm_robust(article_rating ~ article_condition + as.factor(article), clusters = ResponseId, se_type = "stata", data = subset(long_data, long_data$vote_none == 1))

# table

texreg(list(mod1, mod2, mod3, mod4),
       type = "latex",
       file = "results/Table_bias_perception_bypartisanship_vote_articleFEs.tex",
       digits = 3,
       stars = c(0.001, 0.01, 0.05, 0.1),
       symbol = "\\sym{+}",
       include.ci = FALSE,
       omit.coef = "(Intercept)",
       custom.header = list("Article Rating" = 1:4),
       custom.model.names = c("(1)", "(2)", "(3)", "(4)"),
       custom.coef.map = list("article_conditionCNA" = "Central News Agency", 
                              "article_conditionCT" = "China Times",
                              "article_conditionLT" = "Liberty Times", 
                              "article_conditionUD" = "United Daily"),
       custom.gof.rows = list("Respondents" = c("DPP supporters", "KMT supporters", "TPP supporters", "Moderates"),
                              "Article FEs" = c("Yes", "Yes", "Yes", "Yes"), "Individual FEs" = c("No", "No", "No", "No")),
       custom.note = "")

# figure

tidy_mod1 <- tidy(mod1, conf.int = TRUE) %>% mutate(respondent = "DPP supporters")
tidy_mod2 <- tidy(mod2, conf.int = TRUE) %>% mutate(respondent = "KMT supporters")
tidy_mod3 <- tidy(mod3, conf.int = TRUE) %>% mutate(respondent = "TPP supporters")
tidy_mod4 <- tidy(mod4, conf.int = TRUE) %>% mutate(respondent = "Moderates")

combined_models <- rbind(tidy_mod1, tidy_mod2, tidy_mod3, tidy_mod4)
combined_models <- subset(combined_models, combined_models$term != "(Intercept)")
combined_models <- subset(combined_models, !grepl("as.factor", term))
combined_models$term <- factor(combined_models$term, levels = c("article_conditionCNA", "article_conditionCT", "article_conditionLT", "article_conditionUD"), labels = c("Central News Agency", "China Times", "Liberty Times", "United Daily"))
combined_models$respondent <- factor(combined_models$respondent, levels = c("DPP supporters", "KMT supporters", "TPP supporters", "Moderates"))
combined_models$conf.low.90 <- combined_models$estimate - (qnorm(0.95)*combined_models$std.error)
combined_models$conf.high.90 <- combined_models$estimate + (qnorm(0.95)*combined_models$std.error)
combined_models$effect <- "effect"

Figure <- ggplot(combined_models, aes(x = estimate, y = effect, color = respondent, shape = respondent)) + 
  facet_wrap("term", nrow = 1) + 
  geom_point(position = position_dodge(width = 0.5), size = 2) + 
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0, position = position_dodge(width = 0.5)) +
  geom_errorbarh(aes(xmin = conf.low.90, xmax = conf.high.90), height = 0, position = position_dodge(width = 0.5), size = 1) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black", size = 0.75, alpha = 1.0) +
  labs(x = "Estimated Treatment Effects", y = NULL, color = "Respondent", shape = "Respondent", title = NULL) +
  theme(legend.position = "bottom", axis.title.x = element_text(margin = margin(t = 5)), axis.text.x = element_blank(), axis.ticks.x = element_blank()) + 
  scale_x_continuous(limits = c(-2, 2), breaks = seq(-2, 2, length.out = 5), labels = label_number(accuracy = 0.01)) + 
  scale_color_manual(values = c("DPP supporters" = "darkgreen", "KMT supporters" = "darkblue", "TPP supporters" = rgb(0, 0.6, 0.6), "Moderates" = "darkgrey")) + 
  scale_shape_manual(values = c("DPP supporters" = 16, "KMT supporters" = 17, "TPP supporters" = 18, "Moderates" = 4)) + 
  guides(shape = guide_legend(nrow = 1), color = guide_legend(nrow = 1)) + 
  coord_flip()
Figure 

ggsave("results/Figure_bias_perception_bypartisanship_vote_articleFEs.pdf", plot = Figure, width = 10, height = 4.5, units = "in", device = cairo_pdf)
